Abstract¶
Evolve a two-components system and compare the results with the equation (7) in http://link.aps.org/abstract/PRA/v59/pR31. This notebook is meant to be a tool to study pseudo-spin oscillation.
Set up¶
In [34]:
%matplotlib inline
from __future__ import division # to carry out normal division with just one '/'.
import numpy as np
import trottersuzuki as ts
from matplotlib import pyplot as plt
import plotly.plotly as py
import plotly.graph_objs as go
from plotly.offline import download_plotlyjs, init_notebook_mode, iplot
init_notebook_mode()
dim = 100 # square lattice
length = 25.
delta_x = delta_y = length / dim
periods = [0, 0]
kernel_type = "cpu" #the evolution is performed by the CPU kernel
# use harmonic oscillator units
particle_mass_A = particle_mass_B = 1.
w_x = w_y = 1. # frequencies of the harmonic confinement
rot_coord_x = rot_coord_y = dim / 2 # the axis of rotation will be located at the lattice center
def harmonic_pot(x,y):
return 0.5 * particle_mass_A * (w_x**2 * x**2 + w_y**2 * y**2)
TF_radius = 5.938 # wished TF radius, in harmonic oscillator units
coupling_const = np.pi * (TF_radius**4) / 4. # this will be multiplying a wavefunction normalized to unity
# vectors containing the physical coordinates (in h.o. units, and centered in the lattice)
x_vec = y_vec = (np.arange(dim) - dim*0.5) * delta_x
# define the external potential
external_pot_A = np.zeros((dim, dim))
for i in range(0, dim):
y = y_vec[i]
for j in range(0, dim):
x = x_vec[j]
external_pot_A[i,j] = harmonic_pot(x,y)
external_pot_B = external_pot_A
# plot energy and norm
def plot_evolution():
plt.plot(Time, Tot_energy, label = 'Total Energy')
plt.legend()
plt.xlabel('Time')
plt.show()
plt.plot(Time, Norm_A, label = 'Squared norm A')
plt.plot(Time, Norm_B, label = 'Squared norm A')
plt.plot(Time, Norm_A + Norm_B, label = 'Total squared norm')
plt.legend()
plt.xlabel('Time')
plt.axis([0, max_it * iterations * delta_t, 0, np.amax(Norm_A + Norm_B) + 1])
plt.show()
# exports to file the energies at a given time
def get_info(time):
tot_energy = ts.calculate_total_energy_2GPE(p_real_A, p_imag_A, p_real_B, p_imag_B, particle_mass_A, particle_mass_B,
[g_A, g_B, g_AB, Omega_Rabi_real, Omega_Rabi_imag], external_pot_A, external_pot_B,
omega, rot_coord_x, rot_coord_y, delta_x, delta_y)
norm_A = ts.calculate_norm2(p_real_A, p_imag_A, delta_x, delta_y)
norm_B = ts.calculate_norm2(p_real_B, p_imag_B, delta_x, delta_y)
phase_A = ts.get_wave_function_phase(p_real_A, p_imag_A)
phase_B = ts.get_wave_function_phase(p_real_B, p_imag_B)
return [tot_energy, norm_A, norm_B, phase_A[dim/2, dim/2], phase_B[dim/2, dim/2]]
# this is the basic evolution step
def evolution_step():
ts.evolve_2GPE(p_real_A, p_imag_A, p_real_B, p_imag_B, particle_mass_A, particle_mass_B, external_pot_A, external_pot_B,
omega, rot_coord_x, rot_coord_y, delta_x, delta_y, delta_t, iterations,
[g_A, g_B, g_AB, Omega_Rabi_real, Omega_Rabi_imag], kernel_type,
periods, imag_time)
def calculate_mean_ext_pot():
sum_norm_A =0.
sum_norm_B =0.
sum_pot_A = 0.
sum_pot_B = 0.
for yi in range(0, dim):
for xi in range(0, dim):
sum_norm_A += p_real_A[xi,yi]*p_real_A[xi,yi]+p_imag_A[xi,yi]*p_imag_A[xi,yi]
sum_norm_B += p_real_B[xi,yi]*p_real_B[xi,yi]+p_imag_B[xi,yi]*p_imag_B[xi,yi]
sum_pot_A += (p_real_A[xi,yi]*p_real_A[xi,yi]+p_imag_A[xi,yi]*p_imag_A[xi,yi])*external_pot_A[xi,yi]
sum_pot_B += (p_real_B[xi,yi]*p_real_B[xi,yi]+p_imag_B[xi,yi]*p_imag_B[xi,yi])*external_pot_B[xi,yi]
return sum_pot_B / sum_norm_B - sum_pot_A / sum_norm_A
Initialize the state¶
In [35]:
# width of the gaussian envelope
width = 4
def gauss_state(x,y):
val = np.exp(-(x**2 + y**2)/(2. * width**2))
return val
p_real_A = np.zeros((dim,dim))
p_imag_A = np.zeros((dim,dim))
p_real_B = np.zeros((dim,dim))
p_imag_B = np.zeros((dim,dim))
for i in range(0, dim):
y = y_vec[i]
for j in range(0, dim):
x = x_vec[j]
p_real_A[i, j] = np.real(gauss_state(x,y))
p_imag_A[i, j] = np.imag(gauss_state(x,y))
p_real_B[i, j] = np.real(gauss_state(x,y))
p_imag_B[i, j] = np.imag(gauss_state(x,y))
# normalize the initial condition
norm=np.sqrt(delta_x * delta_y * np.sum(np.abs(p_real_A)**2 + np.abs(p_imag_A)**2))
p_real_A=p_real_A/norm
p_imag_A=p_imag_A/norm
norm=np.sqrt(delta_x * delta_y * np.sum(np.abs(p_real_B)**2 + np.abs(p_imag_B)**2))
p_real_B=p_real_B/norm
p_imag_B=p_imag_B/norm
Find the ground state¶
In [36]:
omega = 0.
imag_time = True # True: imaginary time evolution; False: real time evolution
delta_t = 0.5e-4 #evolution time of a single iteration step
iterations = 5000 #number of iterations
start_it=1
max_it = 5
g_A =coupling_const
g_B = coupling_const
g_AB = 0.99 * g_A
Omega_Rabi_real = 0
Omega_Rabi_imag = 0
# plot_evolution variables
Time = np.zeros(max_it - start_it + 2)
Tot_energy = np.zeros(max_it - start_it + 2)
Norm_A = np.zeros(max_it - start_it + 2)
Norm_B = np.zeros(max_it - start_it + 2)
time = 0.
[tot_energy, norm_A, norm_B, a, b] = get_info(time)
Time[0] = time
Tot_energy[0] = tot_energy
Norm_A[0] = norm_A
Norm_B[0] = norm_B
for cont in range(start_it-1, max_it):
evolution_step()
time=(cont + 1) * iterations * delta_t
[tot_energy, norm_A, norm_B, a, b] = get_info(time)
Time[cont+1] = time
Tot_energy[cont+1] = tot_energy
Norm_A[cont+1] = norm_A
Norm_B[cont+1] = norm_B
# calculate overlap between p_A and p_B (only for real p_A and p_B); this is used in the post-processing
k = 0.
for i in range(0, dim):
for j in range(0, dim):
k += p_real_A[i,j] * p_real_B[i,j]
k *= delta_x * delta_y
plot_evolution()
Real time evolution¶
In [37]:
imag_time = False # True: imaginary time evolution; False: real time evolution
iterations = 200 #number of iterations
start_it=1
max_it = 1600
# set Rabi coupling
Rabi_period = 20.
Omega = 2. * np.pi / Rabi_period
delta = np.pi / 3.
Omega_Rabi_imag = 2. * np.pi/Rabi_period * np.sin(delta)
Omega_Rabi_real= 2. * np.pi/Rabi_period * np.cos(delta)
# vectors to be post-processed
Time = np.zeros(max_it - start_it + 2)
Tot_energy = np.zeros(max_it - start_it + 2)
Norm_A = np.zeros(max_it - start_it + 2)
Norm_B = np.zeros(max_it - start_it + 2)
Phase_A = np.zeros(max_it - start_it + 2)
Phase_B = np.zeros(max_it - start_it + 2)
mu_diff_pot = np.zeros(max_it - start_it + 2)
time = 0.
[Tot_energy[0], Norm_A[0], Norm_B[0], Phase_A[0], Phase_B[0]] = get_info(time)
Time[0] = time
mu_diff_pot[0] = calculate_mean_ext_pot()
for cont in range(start_it-1, max_it):
if cont%100 == 0:
print cont
evolution_step()
time=(cont + 1) * iterations * delta_t
[Tot_energy[cont+1], Norm_A[cont+1], Norm_B[cont+1], Phase_A[cont+1], Phase_B[cont+1]] = get_info(time)
Time[cont+1] = time
mu_diff_pot[cont+1] = calculate_mean_ext_pot()
plot_evolution()
Post processing¶
In [38]:
print mu_diff_pot
Calculate $\eta = \frac{N_2 - N_1}{N}$, $\phi_{21} = \phi_2 - \phi_1$, $\frac{d}{dt} \eta$ and $\frac{d}{dt} \phi_{21}$
In [49]:
eta = np.zeros(max_it - start_it + 2)
dt_eta = np.zeros(max_it - start_it + 2)
phase_21 = np.zeros(max_it - start_it + 2)
dt_phase_21 = np.zeros(max_it - start_it + 2)
eta = 0.5 * (Norm_B - Norm_A)
phase_21 = Phase_B - Phase_A
# calculate time derivative of eta and phi_21
for i in range(1,max_it - start_it + 1):
dt_eta[i] = (eta[i+1] - eta[i-1]) / (2 * delta_t * iterations)
var = phase_21[i+1] - phase_21[i-1]
if np.sign(var) * (var) > np.pi:
var = var - np.sign(var)*2*np.pi
dt_phase_21[i] = (var) / (2 * delta_t * iterations)
dt_eta[0]=dt_eta[1]
mu_1 = g_A / length**2 + g_AB / length**2
mu_2 = g_B / length**2 + g_AB / length**2
_dt_eta = -2. * Omega * k * np.sqrt(1. - eta**2) * np.sin(delta + phase_21) # derivative of eta deduced from eta and phi_21
_dt_phase_21 = - (g_A-g_AB) / length**2 * eta + 2 * Omega * k * eta / np.sqrt(1. - eta**2) * np.cos(delta + phase_21)
Plot $\frac{d}{dt}\eta$ (with plotly)¶
In [50]:
x = np.arange(0, 1001 * delta_t * iterations, delta_t * iterations)
trace1 = go.Scatter(
x = x,
y = dt_eta,
mode = 'lines',
name = 'eta from simulation'
)
trace0 = go.Scatter(
x = x,
y = _dt_eta,
mode = 'lines',
name = 'eta from equation'
)
data = [trace0, trace1]
layout = dict(
font = go.Font(size=14),
xaxis=dict(
title='Time',
showticklabels=True,
autotick=False,
ticks='outside',
tick0=0,
dtick=0.5,
ticklen=10,
#tickwidth=4,
),
yaxis = dict(
title = '',
autotick=False,
ticks='outside',
tick0=0,
dtick=0.25,
ticklen=10,
#tickwidth=4,
),
)
fig = dict(data=data, layout=layout)
iplot(fig)
#py.image.save_as(fig, 'eta-time-derivative-null-deltapi3.pdf')
Plot $\frac{d}{dt}\phi_{21}$ (with plotly)¶
In [51]:
x = np.arange(delta_t * iterations, max_it * delta_t * iterations, delta_t * iterations)
trace1 = go.Scatter(
x = x,
y = dt_phase_21,
mode = 'lines',
name = 'phi_21 from simulation'
)
trace0 = go.Scatter(
x = x,
y = _dt_phase_21,
mode = 'lines',
name = 'phi_21 form equation'
)
data = [trace0, trace1]
layout = dict(
font = go.Font(size=14),
xaxis=dict(
title='Time',
showticklabels=True,
autotick=False,
ticks='outside',
tick0=0,
dtick=0.5,
ticklen=10,
#tickwidth=4,
),
yaxis = dict(
title = '',
autotick=False,
ticks='outside',
tick0=0,
dtick=0.25,
ticklen=10,
#tickwidth=4,
),
)
fig = dict(data=data, layout=layout)
iplot(fig)
#py.image.save_as(fig, 'phi_21-time-derivative-null-deltapi3.pdf')
In [52]:
print - (g_A-g_AB)/(length*length)*eta
In [78]:
print Norm_A[250]
In [88]:
x = np.arange(0, 1001 * delta_t * iterations, delta_t * iterations) * 0.1
trace0 = go.Scatter(
x = x,
y = dt_eta,
mode = 'lines+markers',
marker = dict(maxdisplayed = 62),
name = '$ \eta $'
)
trace1 = go.Scatter(
x = x,
y = _dt_eta,
mode = 'lines',
name = '$ f_{\dot{\eta}}(\eta, \phi) $'
)
trace2 = go.Scatter(
x = x,
y = dt_phase_21,
mode = 'lines+markers',
marker = dict(maxdisplayed = 62),
name = '$ \phi $',
yaxis='y2'
)
trace3 = go.Scatter(
x = x,
y = _dt_phase_21,
mode = 'lines',
name = '$ f_{\dot{\phi}}(\eta, \phi) $',
yaxis='y2'
)
data = [trace0, trace1, trace2, trace3]
layout = dict(
legend=dict(
x=0.01,
y=1
),
width = 900,
height = 500,
margin = dict(t=0),
font = go.Font(size=15),
xaxis=dict(
title='$t/T$',
showticklabels=True,
autotick=False,
ticks='outside',
tick0=0,
dtick=0.25,
ticklen=10,
#tickwidth=4,
),
yaxis = dict(
title = '$ \eta $',
autotick=False,
ticks='outside',
tick0=0,
dtick=0.25,
ticklen=10,
#tickwidth=4,
),
yaxis2 = dict(
title = '$ \phi $',
autotick=False,
ticks='outside',
tick0=0,
dtick=0.25,
ticklen=10,
overlaying='y',
side='right',
),
)
fig = dict(data=data, layout=layout)
iplot(fig)
#py.image.save_as(fig, 'phi.pdf')
In [ ]: